﻿using System;
using System.Collections.Generic;
using System.Windows.Markup;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;
using Petra.Gis.ReferenceSystemServices.Transformations;

namespace testSRST
{
  public class testEqTransformJtsk
  {



    public static void testMe()
    {

      //wgs bod v CR
      Coordinate sour = Coordinate.Parse("48°46′50.20133″N, 14°0′49.63044″E");

      //uprava meridianu na ferro a muzeme predstirat ze je to bod na besselove elipsoidu
      //31,6804529;48,7806114805556
      sour.lon -= -17.66666666666666666666666666666667;

      sour = (Coordinate)bessel2Jtsk(sour);

      sour = (Coordinate)jtsk2Bessel(sour);

      //jezek:792355,961974300; 1182894,67248831;
      //lcc:  791835,347664496;-1182117,45616915
      //WGS:  792277.450000000; 1182833.62

    }



    static double alfa = 1.000597498372;
    static double u0 = 49.4599572917;
    static double fi0 = 49.5;

    static double Uk = 59.7118602500 * Math.PI / 180;
    static double Vk = 42.5253936806 * Math.PI / 180;

    static double n = 0.97992470462;
    static double ro0 = 1298039.00462;

    static double kuzpolX;
    static double kuzpolY;

    static double kouleX;
    static double kouleY;

    static double kartX;
    static double kartY;




    public static ICoordinate jtsk2Bessel(ICoordinate coord)
    {
      //sub XY_Bessel(point& P)
      //XY => ro, eps
      //point kart, koule, kuzpol
      //JTSK=P
      kuzpolX = Math.Sqrt(coord.X * coord.X + coord.Y * coord.Y); //ro
      kuzpolY = Math.Atan(coord.Y / coord.X); //eps
      // ro, eps => sirku, delku
      //double n = 0.97992470462, ro0 = 1298039.00462;
      double v = Math.Tan((78.5 / 2 + 45) * Math.PI / 180), t = ro0 / kuzpolX, u = Math.Pow(t, (1 / n));
      kartX = 2 * (Math.Atan(v * u) - (45 * Math.PI / 180)); //sirka
      kartY = kuzpolY / n; //delka
      // sirka, delka => U,V
      //double Uk = 59.7118602500 * PI / 180;
      //double Vk = 42.5253936806 * PI / 180;
      kouleX = Math.Asin(Math.Sin(kartX) * Math.Sin(Uk) - Math.Cos(kartX) * Math.Cos(Uk) * Math.Cos(kartY));//U
      kouleY = Vk - Math.Asin(Math.Cos(kartX) * Math.Sin(kartY) / Math.Cos(kouleX)); //V
      // U,V=> fi, lambda
      //double alfa = 1.000597498372, u0=49.4599572917, fi0=49.5;
      double du = (kouleX * 180 / Math.PI) - u0;

      //bessel lonlat
      coord.X = fi0 + 1.001416022789 * du - 86.87150417 / 1000000.0 * (du * du) +
                16.70197 / 100000000.0 * Math.Pow(du, 3) + 117.5089 / 10000000000.0 * Math.Pow(du, 4); //fi
      coord.Y = ((kouleY * 180 / Math.PI) / alfa); //lambda


      return coord;
    }






    public static ICoordinate bessel2Jtsk(ICoordinate lonlat)
    {
      //point kart, koule,kuzpol,Bessel,JTSK
      //Bessel=coord

      //temporary convert to degrees, because all equations below are in degrees ;-(
      //lonlat.lon = ProjMath.Radians2Degrees(lonlat.X);
      //lonlat.lat = ProjMath.Radians2Degrees(lonlat.Y);

      // Fi, lambda => U, V
      //double alfa = 1.000597498372, u0 = 49.4599572917, fi0 = 49.5;
       double dfy = lonlat.lat - fi0;
      kouleX = (u0 + dfy * (99.8585979496 / (100)) + 86.50351075 / (1000000.0) * (dfy * dfy) -
      15.1091 / (100000000.0) * Math.Pow(dfy, 3) - 117.3673 / (10000000000.0) * Math.Pow(dfy, 4)) * Math.PI / 180; //U
      kouleY = Math.PI / 180 * (alfa * lonlat.lon); //V
      // U,V => sirka, delka
      //double Uk= 59.7118602500 * PI / 180;
      //double Vk= 42.5253936806 * PI / 180;
      //Sirka:
      kartX = Math.Asin(Math.Sin(Uk) * Math.Sin(kouleX) + Math.Cos(Uk) * Math.Cos(kouleX) * Math.Cos(Vk - kouleY));
      //Delka:
      kartY = Math.Asin(Math.Sin(Vk - kouleY) * Math.Cos(kouleX) / Math.Cos(kartX));
      // sirku, delku => ro, eps
      //double n = 0.97992470462, ro0 = 1298039.00462;
      //ro:
      kuzpolX = ro0 * (Math.Pow(Math.Tan((39.25 + 45) * Math.PI / 180) / (Math.Tan(kartX / 2 + 45 * Math.PI / 180)), n));
      kuzpolY = n * kartY; //eps
      // ro, eps => XY
      lonlat.X = kuzpolX * Math.Cos(kuzpolY); //X
      lonlat.Y = kuzpolX * Math.Sin(kuzpolY); //Y

      return lonlat;

    }






  }
}
